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1. INTRODUCTION 

This report treats two important problems in the area of control systems 
design and analysis. The first problem we treat here is the robust stability using 
characteristic polynomial. This problem is treated first in characteristic polynomial 
coefficient space with respect to perturbations in the coefficients of the characteristic 
polynomial, and then for a control system containing perturbed parameters in the 
transfer function description of the plant. In coefficient space, a simple expression is 
first given for the l 2 - stability margin for both monic and non-monic cases. Following 
this, a method is extended to reveal much larger stability region. 

This result has been extended to the parameter space so that one can determine 
the stability margin, in terms of ranges of parameter variations, of the closed loop 
system when the nominal stabilizing controller is given. This stability margin can 
be enlarged by a choice of better stabilizing controller. 

The second problem this report describes is the lower order stabilization 
problem. The motivation of the problem is as follows. Even though the wide 
range of stabilizing controller design methodologies are available in both the state 
space and the transfer function domains, all of these methods produce unnecessarily 
high order controllers. 

In practice, the stabilization is only one of many requirements to be satisfied. 
Therefore, if the order of a stabilizing controller is excessively high, one can normally 
expect to have a even higher order controller upon the completion of design such as 
inclusion of dynamic response requirements, etc. Therefore, it is reasonable to have 
a lowest possible order stabilizing controller first and then adjust the controller to 
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meet additional requirements. 

In this report, the algorithm of designing a lower order stabilizing controller 
is given. The algorithm is not necessarily produce the minimum order controller, 
however the algorithm is theoretically logical and some simulation results show that 
the algorithm works in general. 

The above two problems have been solved and published. These are found 
in Appendix A and B. Finally, some remarks and on going research are briefly 
discussed. 
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2. PUBLICATIONS SUPPORTED BY THE GRANT 
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3. DISCUSSIONS AND DIRECTIONS OF RESEARCH 

The stability margin in the parameter space is in general more practical than 
one in the coefficient space. When parameters enter into coefficients of characteristic 
polynomial linearly, the stability margin is exact. However, when parameters enter 
in nonlinear fashion which is the general facts in practical systems, the margin seems 
to be conservative. 

The results show that the transfer function approach (or polynomial frame- 
work) has some difficulties because in most cases, the state space representation 
gives better description of the dynamic systems in terms of parameters. It is true in 
most cases because the state space description is obtained directly from the dynamic 
equations. 

The next phase of research will be concentrated in developing new idea of robust 
control in the domain of state space. Even though some results of this problem are 
available, they are still primitive and not yet satisfactory. 
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ABSTRACT 

This paper treats the robust stability issue using the characteristic polynomial, 
for two different cases. First in coefficient space with respect to perturbations 
in the coefficients of the characteristic polynomial, and then for a control system 
containing perturbed parameters in the transfer function description of the plant. 
In coefficient space, a simple expression is first given for the l 2 -stability margin 
for both the monic and non-monic cases. Following this, a method is extended to 
reveal much larger stability regions. In parameter space we consider all single input 
(multi output) or single output (multi input) systems with a fixed controller and a 
plant described by a set of transfer functions which are ratios of polynomials with 
variable coefficients. The paper gives a procedure to calculate the radius of the 
largest stability ball in the space of these variable parameters. This calculation is 
important as it serves as a stability margin for the control system. The method is 
based on the application of the orthogonal projection theorem in the appropriate 
Euclidean vector space. The formulas that result are quasi closed form expressions 
for the stability margin, are computationally efficient and provide some insight. The 
paper is illustrated by numerical examples. 
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I. INTRODUCTION 

Kharitonov’s theorem was revealed to the control community in 1984 [1]. 
Since then there has been a growing interest in the theory of robust control with 
structured perturbations, [2]- [6]. Soh, Berger and Dabke introduced the largest 
stability hypersphere around a stable polynomial [2], and in [3] Biernacki, Hwang 
and Bhattacharyya extended their idea to the control problem and introduced the 
largest stability hypersphere in parameter space. This last result is applicable 
whenever the characteristic equation of the closed loop system is a linear or affine 
function of the transfer function subject to perturbation. 

In particular, the latter situation always occurs in single input (multioutput) 
or single output (multiinput) plants, when the parameter vector is chosen to be the 
transfer function coefficients of the plant. These two results, however, suffer from 
the fact that the method given for the calculation of the radius of the stability 
hypersphere is not very satisfactory. In this paper we present new methods for 
this calculation, both in coefficient space and in parameter space, by using the 
orthogonal projection theorem in the appropriate Euclidean vector space. The 
formulas that result are quasi closed form expressions for the radius of the largest 
stability hypersphere. We also present some important simplifications of the scalar 
minimization problem that results. These results constitute a conceptual and 
computational improvement over the stability margin calculations given in [2] and 
[3]. 

The paper also considers the problem of computing the l°° margin around 
a Hurwitz polynomial and and interesting method is given which is based on 
Kharitonov’s theorem and which greatly simplifies the one given in [1] and then 
[5] . Moreover it is shown by an example that this procedure is very simply extened 
and yields stability regions which are much larger than those provided by the simple 
computation of the / 00 -margin. The paper is organized as follows. In section II, we 
present the results concerning the computation of stability margins in coefficient 
space, and in section III we introduce the general setting for the control problem 
considered in this paper, and present our result for the calculation of the l 2 -stability 
margin. 


II. STABILITY MARGINS IN COEFFICIENT SPACE 

In this section we first consider a stable polynomial S(-) of degree n, and we 
give an expression for the radius, in the space of coefficients, of the largest stability 
hypersphere around the nominal point, both in the non-monic and monic case. 


II. 1. I 2 stability margin 
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a) The non-monic case: 

Let S(-) be an arbitrary stable polynomial of order n, 

8(s) = + 8\S + ‘ ' * + 8 n S n . 

The l 2 norm of £(s) is defined by, 




1=0 


Let us now separate 8(s) into its odd and even parts: 


8(s) = 

even degree terms odd degree terms 

and let us also define 8 e (to) and 8° (to) as follows: 


8 e ( w) = 8 even (ju>) = 8 0 - 8 2 u 2 + 8 4 u> a - 
8° (u>) = =6 X - 8 3 to 2 + 8 5 to 4 - 


( 2 . 1 ) 


Theorem 2.1 

The radius of the largest stability hypersphere around £(s) is given by: 


p(8) = min(|£o|, |£»l, inf ’ d%), 


where d" is given by: 
i) n = 2p 

d n2 _ 

w 1 + u; 4 H +w 4 P i "l+o; 4 H fu; 4 ^" 1 ) 

ii) n = 2p + 1 

,„2 = [£Ml+i£Ml ! 

“ l + w 4 h 


( 2 . 2 ) 


( 2 . 3 ) 


b) The monic case: 
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In this section, we are given a stable monic polynomial 
0(s) = Po + ft\S + + 0n$ n 1 + 

and we want to find the radius of the largest stability hypersphere in the affine 
space of minoc polynomials of degree n. We have the following result. 

Theorem 2.2 


p{P) = min(|/? 0 |, inf d"') 

oj>0 


where 

computed for 


that is 


d 


nf 


= d 


n 

u 


f3(s) - s n - u 2 s n ~ 2 


i) n = 2p 


d „,2 = [<5‘ M ] 2 + |5»] 2 

" 1 +u> 4 + ... 


ii) n = 2p + 1 

d n / 2 _ [^ M ] 2 [^ V )] 2 

1 + CO 4 + • • • + U> 4p 1 + W 4 + • • • + u>4(p-i) 


(2.4) 


(2.5) 


we now show how the minimization problem that results from the application of 
formulas (2.2) — (2.5) can be simplified. Consider for example the calculation of 


dmin — inf d w . 

w>0 


A simple manipulation we show that there is no need to carry out a minimization 
over the finite range [0, oo). We will consider the case when n = 2 p, but a similar 
derivation holds if n is odd. 

First it is clear that 


dm in=min( inf d 2 u , inf d\), 

wG[0,l] w€[0,l] ^ 
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and then we have 

l<5»i 2 = (s 0 - s 2 u? + • • • + (-1 ) p s 2 yn 2 

|6°(u;)| 2 = (^i — 830!" + • • • + (— l) p 1 62p-iu> 2p 2 ) 2 , 

which yields: 

= ^(hp - Sir-iJ + ■ ■ • + 

and 

l^°(~)| 2 = a ,4(p-I) (^2p-l - t>2p- 3 ^ 2 H h (-I/ -1 t>\U 2p 2 ) 2 . 

So that finally, 

2 _ ■£?(& 2p ~ ^2p-2^ 2 -I 1" ( l) P ^QUJ 2p ) 2 

l + ^T + ^+-” + ^F 

^(p-iy C^p-i -^ 2 p- 3 ^ 2 H 1- (-l) p ~ 1 ^io> 2p ~ 2 ) 2 

1 + "* +* u.-'fp- 1 ) 

This last expression however is nothing but 

^2 _ (^2p ~ ^2p-2<* >2 + <^ 2 p- 4 ^ )4 + • • • + ( — 1) p <$qO> 2p ) 2 

i 1 -f a; 4 + u> 8 + w 4p 

(<^2p-l ~ <^2p— 3 <^ 2 + <^ 2 p- 5 ^ 4 + • ‘ • + ( — l) ;) ~ 1 l 5 ik; 2p ~ 2 ) 2 

1 + u> 4 + u> 8 H + u> 4 (p-i) 

and we can see that 8\_ has exactly the same structure as d 2 . Can S 2 ^ be considered 

(jJ u> 

as the “<£” of some other polynomial? The answer is of course yes. Consider 
8'(s) = s n d(j), which in our case is 

S'(s) = s 2p 8( — ) = 6 2 p + <!> 2 p_is + • • * + 8qs 2j> , 
s 

then clearly 

d ,e (u>) = 6 2 p — ^2p-2^ 2 + • • • + ( — l) P 8()U 2p 

and 

6' (u>) = 8 2p —i — <^2p— 3U; 2 + • • • + ( — l) p 1 8\uj 2p 2 . 
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Thus we see that in fact 8\ corresponds to d £ computed for d'(-). Suppose now that 

UJ 

you have a subroutine DM IN (8) that takes the vector of coefficients d as input and 
returns the minimum of 8^ over [0, 1]. Then the following algorithm will compute 
d m in by simply calling DM IN two times. 

1 . Set d ( 8q , 8 \ , • * ■ , 8ji ) . 

2. First call: d x — DMIN(8). 

3. Switch: set 8 = (d„, d n _i , • • • , 8 0 ). 

4. Second call: d 2 = DMIN(8). 

5. d min = mm(d x ,d 2 ). 

Incidentally, we already knew that the two polynomials 8(s) and 8'(s ) = s n d(j) 
are stable together (i.e. one is stable if and only if the other one is stable). The 
development above tells us that moreover p(8) = p{8'). 

It is to be noted that the results of this section were first proved in [4] . Recently 
similar results have been achieved in [8]. We now turn to the problem of finding 
the /°°-stability margin. 

II. 2 l°° stability margin 

We consider here as given a stable nominal polynomial d(.s), as well as a set 
of non-negative numbers a 0 , a\, • • • , a n , and we want to find the largest box of the 
form 


B p = (do - &oP, do + Oiop) x (di - a x p, 8 X + at\p) x • • • x (d n - a n p, d„ + a n p ). (2.6) 

Kharitonov’s theorem tells us that the closed box, 

B'p = [d 0 — c*oP,do + olqp] x [di - di -(- a x p\ x ■ • • x [d n - a n p,8 n + a n p ]. (2.7) 

is stable if and only if its four Kharitonov polynomials are stable. Denoting 

K even (s) = Oi 0 — a 2 s 2 + a 4 s 4 , 

A' odd (s) = ais — a 3 s 3 + , 

and as in (2.1), 

K e (u ) = n 0 + a 2 uj 2 + 0:4a; 4 4 , 

K°(u>) = ai + a 3 u 2 -f a 5 u ; 4 -I , 


(2.9) 
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the four Kharitonov polynomials associated with B' p can ben written as, 

A» = S(s) - pK™*(s) - pK M (s\ 

K%s) = S(s) - pK'™(s) + pK odd (s), ^ 

K%s) = Ht) + pK*™(s) - pK^is), 

A» = *M + pK’™(s) + pI< M (s). 

THus we know that when increasing p, B' p will first contain an unstable polynomial 

or a polynomial of degree < n when one of the four polynomials above becomes of 
degree < n or acquires a ju> root or a root at the origin. 


The case of a root at the origin or the case of a loss in degree is trivially achieved 
for 


P = 


<*o 


and 



respectively. 


( 2 . 11 ) 


Let us now look at the case of the appearance of a ju> root for u> > 0. If we 
consider for example K p (s), we know that this polynomial has a ju> root if and only 
if 

r = s’M - pK e (u>) = o 

{ and (2.12) 

( if*» = «”(“)- pK‘(u) = 0. 

This of course is possible if and only if we have 

8 e {u)I<°(u) - 8°(uj)K e (uj) = 0. (2.13) 

THis same polynomial is associated with K p (s ) which has a ju> root if and only if 

r K 4 p ‘(u;) = 6 e (u;)+pK e (u) = 0 

< and (2-14) 

i x-‘» = «>) + a>a'» = o. 

Note that this polynomial is a polynomial of degree n - 1 in w 2 and therefore we 
just have to find the positive roots of a polynomial of degree n — 1. Once these 
roots are found, we then look at the values (at these points) of the ratios 

>(u;) 6°(u>) 

~ K°(lo) ' 



s 


The minimum positive value of these ratios is denoted p\ and corresponds to K 1 (s), 
and the negative of the maximum negative value is denoted p 4 and corresponds to 
Iv 4 (s). 

Similarly, the polynomial in u > 2 associated with A' 2 (s) and A' 3 (s) is 

S e (uj)K 0 ( w) + S°(u)K e (u) = 0. (2.15) 

and p 2 , P 3 are defined in a similar fashion. The maximum centered box is then 
obtained by taking p to be equal to the smallest the six positive numbers 

and 

<* 0 «n 


Once this first box is obtained it is also possible to further extend it, if the 
center is allowed to move, by freezing the coefficients corresponding to the ‘closest’ 
I\'(s ) in the first stage and using the same approach. An example will best explain 
the details of this extension procedure. 

Example : 


Consider the case treated in [7] of the polynomial, 

S(s) = s 6 + 14.0s 5 + 80.25s 4 + 251.25s 3 + 502.75s 2 + 667.25s + 433.5, 
and the set of parameters, 

<*o = 92.32, a 4 = 33.36, a 2 = 38.28 
c*3 = 15.075, «4 = 6.2, e*5 = 1.4, a 6 = 0.14 

Here we have 

K\u) = 92.32 + 38.28u> 2 + 6.2a> 4 + 0.14u> 6 
K°(u) = 33.36 + 15.075a; 2 + 1.4u; 4 . 

Letting t = o> 2 , we have to find the positive roots of the two following polynomials, 

Pi(t) = 6 e (t)K°(t) - S°(t)K e (t ) 

= -47138.96 - 12583.6575* - 106.49625* 2 + 1400.97375* 3 
+ 45. 65* 4 — 3.36t 5 . 
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and 

P 2 (t)=6 e (t)K°(t) + 6°(t)K e (t) 

= 76062.08 - 7889.7975* - 8483.33625< 2 - 455.856254< 3 
+ 148. 9* 4 + 0.56* 5 . 

It turns out that Pi(t) has two positive roots, 

t{ = 4.020612 and t\ = 28.1620029, 
and P 2 {t) also has two positive roots which are 

t\ = 2.5415548 and t\ = 9.0241863. 

Based on these roots one computes the following values for the pi's, 

Pl = 2.9937539, p 2 = 1.6229978, 

p 3 = 1.4757364, p 4 = 1.0001038. 

From this we conclude that the l°° - margin is obtained by putting p = p 4 = 
1.0001038 in (2.7) giving, 

B[ = [341.17042,525.82957] x [633.88654,700.61346] 
x [464.46603,541.03397] x [236.17344,266.32656] 
x [74.04936,86.45064] x [12.59986,15.40014] 
x [0.85999,1.14001]. 

We know that the Kharitonov polynomial K 4 (s) associated with B\ is ‘unstable’. 
However we can still increase B[ in the other direction and thus consider, 

B[ p , = [£ 0 - otop 1 , 525.82598] x [S 4 - a lP ', 700.61346] 
x [464.46603, S 2 + a 2 p'\ x [236.17344, S 3 + oc 3 p'] 
x [S 4 - a 4 86.45064] x [S 5 - a 5 p', 15.40014] 
x [0.85999, <$6 + a&p']- 

The 3 remaining Kharitonov polynomials associated with B' lp , can be rewritten as 

I<l,(s) = 6(s) - p'K e (s) - p'K’(s), 

K}(3) = 6(3 ) - p'K e (s) + pK°(s), 

Kp(s) = 6(s) + pK‘(s) - p'K°(3). 
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The smallest p' for which I\ pl (s) gets a joj root is of course the same as the one 
calculated above, that is p\ = 2.9937539. 

For Kp>(s ) on the other hand, p' 2 is simply given by the smallest positive value 
of the ratio 6 e (u)/ K e (uj), at the real roots of, 

S°(u) + pK°(u>) = K 4 o (uj) 

= 700.61346 - 236.17344a; 2 + 15.40014a; 4 . 

this gives the value, 

p' 2 = 2.3459473. 

Similarly, p' 3 is given by the smallest positive value of the ration S°(uj)/ K°(uj), at 
the real roots of, 

6 e (u) + pK e (u) = K ie (to) 

= 525.82958 - 464.46603a; 2 + 86.45064a; 4 - 0.85999a/ . 
and this yields the value 

p' 3 - 4.918607463. 

Hence, the value of p' is p' = p 2 = 2.3459473. This gives rise to the augmented box, 

B 2 = [216.92215,525.82958] x [588.98920,700.61346] 
x [464.46603,592.55286] x [236.1734,286.61515] 
x [65.70513,86.45064] x [10.71568,15.40014] 
x [0.85999,1.32843]. 

Here again we know that the Kharitonov polynomials K 4 {s) andA' 2 (s) associated 
with B 2 are ‘unstable’. However we can still increase B 2 by considering, 

B' 2p „ = [216.92215,525.82958] x [<5i - a lP ", 700.61346] 
x [464.46603,592.55286] x [236.17344, S 3 + a 3 p''] 
x [65.70513,86.45064] x [$ 5 - a 5 p", 15.40014] 
x [0.85999,1.32843]. 

The two remaining Kharitonov polynomials are 

K' p „(s) = 6(s) - p'K\s) - p"K°($), 

1 <1 „(i <) = 6(s) + pK‘(s)-p"K°(s). 
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The value of p'' 3 can be seen to the same as p 3 that is 

p 3 = 4.918607463. 

As for p'\ it is given by the smallest positive value of the ratio 6(u))/K°( a;), at the 
real roots of, 

6 e (uj) - p'K e (u) = 

216.92215 - 592.55286a; 2 + 65.70513a; 4 - 1.32843a; 6 . 

This yields p'[ = 4.287652665, which therefore gives the final answer, 

p" = p'\ = 4.287652665, 
so that we get another increase and a final 

B' 2p n = [216.92215,525.82958] x [524.21391,700.6134] 
x [464.46603,592.55286] x [236.17344,315.88636] 
x [65.70513,86.45064] x [7.99729,15.40014] 
x [0.85999,1.32843], 

for which three of the four associated Kharitonov polynomials, namely 7T 1 (s), K 2 (s ) 
and 7f 4 (s) are ‘unstable’. 

We now look at the control problem and consider the stability margin in 
parameter space. 

III. THE / 2 - STABILITY MARGIN IN PARAMETER SPACE 

In the standard feedback system of Figure 1 , suppose that the plant is either 
single input (multioutput) or single output (multiinput). Since the formulation is 
similar for the two cases we restrict our considerations to the single input case. 


tC 


M «■ 






Figure 1 : Unity Feedback System 
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Let therefore, 


G« = 


dP(s ) 


"?(*) 




be the transfer function of the system. The order of the plant is q and 


d p (s) — d p s q + • • • + d p , 

while 

"?(*) = n lq s9 + •'• + <, o- 

The controller transfer function on the other hand is of order r and is described by 

C( s ) = ' ' ' n "*^’ 

where 

d c (s) = d c r s r + --- + d c 0 , 

n i( s ) = n lr sr + -” + < 0 - 

The characteristic polynomial of the closed loop system is then given by the 
polynomial S(-) of degree n = q + r 

S(s) = d c (s)d p (s) + n c m (s)n p m (s ) + • • • + nj(s)n?0). (3.1) 

In this paper the plant parameter vector is taken to be 

P:=[n?(s) n p (s) ■ ■ ■ n^s) <P(s)]. 

The purpose of this section is to derive a measure of stability (stability margin) for 
plants with perturbed coefficients. This can be done by finding the largest stability 
hypersphere in parameter space. We give here a new method for computing the 
radius of such a sphere. For any / and n we denote by V l n the vector space 


V n xV n x---xV n 

' ^ ' 

/ times 

where x designates the Cartesian product and V n is the vector space of real 
polynomials of degree less them or equal to n. V l n is supplied with the natural 
induced inner product defined as follows: 
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If 

£=[AW P 2 (s) ■ • • P,(3)], 

and 

R=[Ri(s) R 2 (s) • • • R,(s)) t 

then 

/ i n 

« P,R »= £ < PM,Ri(s) >= ££p M r,„. < 3 - 2) 

i=l 1=1 j = 0 

The norm associated with this inner product corresponds to the Euclidean norm. 
As an example, for a plant of order q as defined above, we have 

£ = KW «5W • • • <(*) <*' W]sP” +1 , 

and 

m q q 

iiai 2 = ££</ + £<• 

1=1 j=0 jf = 0 

Here again, the question arises of being able to define a stability margin, but this 
time in parameter space. This stability margin will then tell us how much we can 
perturb the original plant P and yet remain stable. The largest stability hypersphere 
as it was defined in [3] or [9] is characterized by the following theorem. 

Theorem 3.1 

Let P = [nj(s) • • n£j(.s) d p (s) ] be a given plant of order q, and C(s) 

a stabilizing controller of order r described by X_ = [nj(s) • • n^s) d c (s)], 

and let n = q + r. 

a) There exists a largest stability hypersphere S(P_,p(P_)) around P, which is 
characterized by 

• For every plant P' within the sphere, the closed loop characteristic 
polynomial (*> X (P') is stable and of degree n. 

• At least one plant P i' on the sphere itself is such that 6 x (Pi!) is unstable 
or of degree less than n. 

b) Moreover if P^ is any plant on the sphere such that 8 X (P 1 ') is unstable, then 
the unstable roots of 8 X (P'') can only be pure imaginary or zero. 

The radius p(P) of the largest stability hypersphere around P is now given by 
the following result. 
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Theorem 3.2 

Let P = [nf(.s) • • • n£j(.s) d p (s)] t> e a plant of order q and let C(s ) be a 

stabilizing controller of order r determined by 

X = [nf(s) ••• n c m (s) d c (s)]. 

The radius of the largest stability hypersphere around P is given by 

p(P) = min( inf d p , d p , d p n ), 

w>0 

where 

a) 

,„ 2 _ A?[[^ 1 ]] 2 + Aj[[Z 2 ]] 2 -2A 1 A 2 «Z 1 ,Z 2 » 

M 2 [£ 2 ]] 2 -<<£ i ,£ 2 >> 2 

with 

m 

Ax = ^(n? e (u;)n?V)-w 2 nrM<V)) 

i=l 

+ d ce (u)d pe (cj) - u 2 d co {uj)d po { w), 

m 

^ = S(nr(«w(‘»)+"“M"r(o’)) 

i=i 

+ d ce (u)d po (u;) + d co (u)d pe (u), 

and 

2i = K e (u;)Pi( S ) + nJ»P2(5), ,d ce ( W )Px(s) + d co (a;)P 2 ( S )), 

Z 2 = (nr(«)ft(«) -« 2 nfV)Ato, , d c »P 2 (s) - a, 2 d c »Px(s)) 

where Pi(-), P 2 (-) are defined as follows (depending on the plant order 5): 

i) q — 21 

p / \ _ / 5 — u; 2 s 2 H h (— l) ,-1 u? 2,-2 s 2i_1 if <7 ^ 0 

^ lW- \0 if 9 = 0 

P 2 (s) = 1 — -| 1- (— 1 ) l u 2 l s 21 . 


ii) q = 2/ + 1 
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b) 


c ) 


Pi (a) = s-o> 2 s 3 + ■■■ + (-l) l u 2l s 2l+ \ 
P 2 (s) = 1 - u 2 s 2 H 1- (— l) / u; 2, .s 2 ^ 


do 


2 




12^=1 n i,o ra ;,o 


V m n<? 5 
Z_a=i 'Vo 


+ «) 2 

+ ^§ 2 


<C = 


(Eli 


77 • 77^ 

n i i q n t 1 r 


+ dw r y 


T m n? : 

Z^i=l '*t,r 


+ 4 


The proof of this result relies on the application of the projection theorem in 
the Euclidean vector space V™ +1 and can be found in [9]. We now provide an 
example to show the applicability of this result. 

Example : 


Consider the following single input, single output plant of order q = 3, 

r(A = nP M = 5 

W dP(s ) 1 - a + 4s 2 + a 3 

A stabilizing controller for G(s) is 

c(s) = nIM = _JL_ 

C{ } d c (s ) 1 + a 

which is of order r = 1 and the resulting characteristic polynomial is 

S(s) = n c (s)n p (s) + d c (s)d p (s) 

= l + 3a + 3a 2 +5a 3 + a 4 . 

In this case we immediately have 

1 


do 2 = jp and d p n 2 = 1. 


And for , we compute 


2 3 


Pi(s) — s — u rs 


2 


P 2 (^) = 1 — ars 


and 


5 
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n 


\u) = 3, n co (u>) = 0, d ce ( u) = 1, d co {u) = 1, 


n 


pe/ 


and then 


and 


so that 


and also 


(w) = 0, n yo (u) = 1, d pe { u) = 1 — 4ar , d p0 (o;) = -1 - u/ 

Ai = <5 e (o;) = 1 - 3a; 2 +w 4 , A 2 = 6°(u) = 3 - 5a; 2 , 

Z 1 =(ZP 1 (s),P 1 (s) + P a (s)) 

= (3s — 3a; 2 s 3 , 1 + s — a; 2 s 2 — a; 2 s 3 ), 

Z 2 = (3P 2 (s),P 2 (s)-a; 2 P 1 (s)) 

= (3 - 3a; 2 s 2 , 1 — a; 2 s - a; 2 s 2 + a; 4 s 3 ), 

[[Z,]] 2 = 11 + 11a; 4 

[[ 2. 2 ]] 2 = 10 + Hk ; 4 + a ; 8 , 


<< Z_ X ,Z_ 2 >>= 1 — a; 2 + a; 4 — a; 6 = (1 — a; 2 )(l + a; 4 ). 


Finally we get 


d p = 

U u; 


11(1 - 3a; 2 + a; 4 ) 2 + (10 +a; 4 )(3 - 5a; 2 ) 2 - 2(1 - a; 2 )(l - 3a; 2 )(l - 3a; 2 + a; 4 )(3 - 5a; 2 ) 

(1 + a; 4 )(109 + 2a; 2 + 10a; 4 ) 

And the minimization of this function over [0, +oo) yields d p 2 ~ .012678. 

Some More Remarks : As before in the coefficient space, it is possible to carry 
out two similar minimizations on the finite range [0, 1] in order to compute 
^min = inf w <o d p . After one minimization, one operates the following permutations. 

n^(s),d p (s) are replaced by s 9 nf(-) and s 9 c? p (-) 

s s 


whereas 

n?(s), d c (s) are replaced by s r n?( — ), s r d c (-). 

s s 

For example in the case that we treated above, one just replaces 
n c (s) = 3 by n' c (s ) = 3s, whereas 
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d' c (s ) = d c (s ) is unchanged 
and 

n p (s) = s by n' p (s) = s 2 and 

d p (s) = 1 — s + 4 s 2 + s 3 by d' p (s ) = 1 + 4s — s 2 4- s 3 . 
Example 2 : 


The Figure shows the experimental model which was developed by NASA 
Langley Research Center in order to demonstrate slewing flexible structures in a 
single axis while simultaneously suppressing vibration motion by the end of the 
maneuver. The detailed information of the model may be found in [10]. If we 
consider that the stiffness coefficient El and motor viscous drag c are uncertain 
parameters, we can redefine parameters as follows. 


d:= ( k t k b /R a +c)N 2 

e := (1.875) V^/(/>0 4 


Now we introduce the zero-th order controller 


C(s) = [raj «2 ri 3 n 4 


Then the characteristic equation becomes 


6(s) = s 4 

+ (0.032394788d+ 1.688950471n 3 + 21.95051139n 4 )s 3 
+ (1.031307514e + 1.688950471^ + 21.95051139n 2 )s 2 
+ (0.032585741de + 1.69806062en 3 )s 
+ 1.69806062en-! 

= 0 



ANGULAR POTENTIOMETER 
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We can simply rewrite the equation to 

S(s) = d (0.032394788s 3 ) 

Rd») Ql ( a) 

+ e (1.031307514s 2 + 1.69806062n 3 s + 1.69806062^) 

' 

Kj(s) Qa(«) 

+ de (0.032585741s) 

' ' 

Qa(») 

+ B(s ) 

= 0 


where 

S(s) = 

s 4 + (1.688950471n 3 + 21.95051139n 4 )s 3 + (1.688950471^ + 21.95051139n 2 )s 2 
Thus 

6(s) = i?i(s)Qi(s) + R,2(s)Q2(s) + R 3 (s)Q 3 (s) + B($) 

The value of d v n is in fact +oo since no matter what controller you choose, 
the characteristic polynomial remains of order 4. The same argument can be made 
for d v 0 . As long as the controller keeps the 0 l h and the highest order coefficients 
of the characteristic equation positive, these two coefficients will not effect the 
stability of the characteristic polynomial. Therefore, d £ determines the stability of 
the polynomial. In order to determine we have 


and, 

Then, 


Pi(*) = o, 


P 2 (s) = 1, 


Xi = S e (ui) 

= 1.6S98906162en x - (1.031307514e + 1.688950471^ + 21.95051139n 2 )a; 2 

+ w 4 


\ 2 =8°{u) 

= (0.032585731de 

+ 1.69806062en 3 ) - (0.032394788d + 1.688950471n 3 + 21.95051139n 4 )u; 2 
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«i = («i(").<K(“).OS(“))> 

Za = (<?!("). 05("). «(“))• 

Thus, 

Z x = (-0.032394788a; 2 , 1.69806062n 3 , 0.032394788), 

Z 2 = (0, 1.69806062m - 1.031307514a; 2 , 0). 

[[Z a ]] 2 = (-0.032394788a; 2 ) 2 + (1.69806062n 3 ) 2 + (0.032394788) 2 , 

[[Z 2 ]] 2 = (1.69806062m - 1.031307514a; 2 ) 2 . 

<< Z\, Z 2 >>= 1.69806062n 3 (1.69806062m - 1.031307514a; 2 ). 

Then, we can use d v w formula. The initial stabilizing controller was chosen to 
be 

C(s) = [0.0369687 -16.647767 -29.677343 0.449643] 
and it provides 

d p = 7.62351163 

UJ 

= allowable perturbation of \/(Ad) 2 + (Ae) 2 + (A(t/e)) 2 

The obtained robust controller is 

C*(s) = [369.73937 113.46374 48.22846 8.208517] 
and it provides 

dl = 10248.294586 

= allowable perturbation of \Z(Ad) 2 + (Ae) 2 + (A(de)) 2 

This result is expected to be somewhat conservative because the nonlinear term 
de was treated as the separate independent parameter. 
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ABSTRACT 

This paper presents an algorithm for stabilizing a linear multivariable system 
with a controller of fixed dynamic order. This is an output feedback stabilization 
problem. An algorithm attempts to solve this via a sequence of approximate 
pole assignment problems. The approximation is driven by the optimization of 
a performance index consisting of a weighted sum of the condition number of the 
closed loop eigenvectors and the norm of the difference between the computed and 
actual controls. 

The algorithm can be used for generating low order solutions to the regula- 
tor problem. The problem treated here is useful in design problems which involve 
parameter optimization and is also important in practical situations where stabi- 
lization is to be accomplished with a fixed number of available parameters. 



2 


1. INTRODUCTION 

The regulator or feedback stabilization problem is the basic problem that 
control theory attempts to solve. Many design procedures can only be initiated 
after a nominal stabilizing controller has been found. However, except for very 
special cases, there axe no direct procedures available to solve this problem when 
the controller order is fixed. Existing solutions to the regulator problem can only 
generate controllers that axe of high enough order that arbitrary pole placement 
becomes possible. This includes the LQG theory [1], observed state feedback [2] 
and arbitrary pole placement approaches [3] [4]. Controllers that are robust with 
respect to unstructured perturbations evidently suffer from the same difficulty of 
high order (see examples given in [5]). We also mention that adaptive control theory 
is notorious for producing high order solutions. 

It is certainly essential in practice, to have low order solutions to the stabi- 
lization problem. This requirement arises because the controller must eventually 
carry out several functions such as tracking, disturbance rejection, desensitization 
against parameter variations, provide good transient response, small steady state 
error, prevent various signals from saturating etc., in addition to the basic task of 
stabilization. Many of these requirements are in conflict with each other in ways 
that cannot be handled analytically and the only recourse left to the designer is 
to iteratively redesign the controller using adhoc methods and graphical displays 
until a satisfactory solution is obtained. This redesign must be carried out in the 
parameter space of the stabilizing controller. If the basic stabilizing controller order 
is unnecessarily high this parameter space is also of high dimension and the subse- 
quent design process can become unwieldy. From this prospective, the high order of 
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controllers produced by “modern” control theory is one of the severest limitations 
of this theory. 

We attempt to alleviate this problem by presenting, in this paper a direct 
algorithm in the state space domain, for designing low order stabilizing controllers. 
This algorithm first attempts to stabilize the closed loop system with a fixed order 
controller. This corresponds to an extended output feedback stabilization problem 
for which no analytical solution is available. We attempt to solve this iteratively. 
At each iteration a state feedback matrix assigning a prescribed set of eigenvalues is 
found and this matrix is approximated by output feedback. This is done successively 
by readjusting the desired closed loop pole locations in the left half of the complex 
plane to minimize a performance index that measures the deviation of the actual 
eigenvalues from the desired ones. A low order solution is found by sequentially 
increasing the controller order until stabilization is achieved. 

The algorithm that is given depends on the parameterization of the state 
feedback pole assignment problem derived in [6]. This is briefly described in the 
next section. In Section 3, the fixed order output feedback stabilization problem is 
formulated as an optimization problem and Section 4 describes how the performance 
index can be decreased by increasing the controller order. Examples are given in 
Section 5 and some of the gradient evaluations of Section 4 are derived in the 
Appendix. 

2. THE SYLVESTER EQUATION FORMULATION 

An algorithm was introduced in [6] for solving the pole assignment problem 
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using state feedback. This algorithm consists of solving for X and the for F 


AX - XA = -BG (2.1) 

FX = G (2.2) 

for given (A, B, A) with an arbitrary choice of G. In (2.1) and (2.2) A , X and A 
are nxn matrices. From a result in [7] the solution X of (2.1) generically has full 
rank if (A, B) is controllable and (G,A) is observable. Let A ,(T) denote the i ih 
eigenvalue of T and A(T) the spectrum or eigenvalue set of T. It follows that if X 
has full rank the solution F has the property: 

A (A + BF) = A .(A) (2.3) 

The advantage of this algorithm are: 

a) The algebraic variety F( A) of matrices F which assign a prescribed set of 
eigenvalues A can be obtained by setting A = A(A) for a fixed A , and letting 
the free parameter G run through the set of all possible real values. 

b) Efficient numerical procedures [8] are available for the solution of Sylvester’s 
equation (2.1). 

Based on this parameterization of F( A) algorithms were given [9] and [10] for 
optimizing the conditioning of the closed loop eigenvectors and [11] for minimizing 
the norm of the state feedback matrix F. Here, we extend these results by 
considering measurement rather than state feedback and by treating the problem 
of stabilization rather than arbitrary pole placement. 


3. OUTPUT FEEDBACK CONTROLLERS 
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Consider the linear time invariant plant S cascaded with the p th order feedback 
compensator C. 

S : x = Ax + Bu 

(3.1) 


The closed loop system is 


x 

x c 


S : x = Ax + Bu 
y m = Cx. 

C : x c = A c x c + B c y m 
u “ C c x c “I - D c y m 


A + BD c C BC c 
B c C A c 


x 

X c 


(3.2) 


(3.3) 



and the transfer function of the p th order compensator is 

C(s) := C c (sl - A c )~ l B c + D c (3.5) 


The formula (3.4) shows that any fixed order compensator design problem 
is equivalent to a static output feedback problem. In particular the problem of 
stabilization with a fixed order controller p is equivalent to that of stabilizing 
A p + B p K p Cp by choice of K p . The general solution of this problem is unknown. 
The best available special results axe those of Brasch and Pearson [3] and Kimura [4] 
which deal respectively with arbitrary eigenvalue assignment and “almost” arbitrary 
eigenvalue assignment. 


Let A denote a symmetric set of n+p complex numbers (i.e. complex numbers 
occur in complex conjugate pairs) and let 

K r ( A) := {Kp\K p € R , "' +r) " r+ ' , Kot A r + B,K P C P ) 6 A) 


(3.6) 
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where A p <E #("+?)*("+?), Bp e R (n+ P )x(m + P ) ? and ^ € ^(r+ P )x(n+ P ) are ag in (3 4) 

The result of Brasch and Pearson [3] states that if (A, B,C) is controllable 
and observable with controllability index v c and observability index i/ 0 , and p > 
min{i/ c , u 0 } then K p (A) ^ 0 for every choice of A. The result of Kimura [4] states 
that ifp>n — m — r + 1 then a{A p + B p K p C p ) can be made arbitrarily close to 
any set A of n + p symmetric complex numbers. 

The lower bound on the order of a stabilizing controller established by the 
above results is in general too conservative. This stems from the fact that both 
results essentially require arbitrary pole placement. In fact for specific choices of A, 
K p ( A) will “almost always” be empty unless p the compensator order is high. To 
lower the compensator order we therefore relax the specification of A in (3.6) to a 
simply connected region ft C C~ and consider the family 

£ p (ft) = {I< p \K p e ij(n*+ P )x(r+ P ) 5 A(Ap + Bp K p C p ) C ft C C~}. (3.7) 

It is reasonable to expect that K_ p {Vt) will in general be nonempty for values of 
p much less than the lower bounds given by the results of Brasch and Pearson or 
Kimura and numerical examples support this. 

The effective characterization of the family K p (ft) is an unsolved open problem. 
Our approach to this problem will be to consider the state feedback family 

Z p (ft) = {Fp\F p E A (Ap + B P F P ) C ft C C~} (3.8) 

and determine an F p € Z p (ft) and then find K p such that ||.F p — K p C p \\ is small 
in the hope that such a K p 6 2£ p (ft). The advantage of this approach is that the 
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family F p (Q.) can be characterized conveniently as shown later. For the remainder 
of this section we drop the subscript p for convenience. 

In general, even if \\F — KC\\ is small it is not in general true that A (A + BF) 
and \(A+BKC) are close. The latter can be achieved by making the eigenstructure 
of A + BF as orthonormal as possible. Let cr moz (T) and cr mtn (T) denote the largest 
and smallest singular values of T. It is well known [8] [12] that the perturbation of 
the eigenvalues of the diagonalisable matrix (A + BF) for changes in the entries is 
small if the condition number k(X) := ||X|| 2 ||X -1 ]| 2 of the eigenvector matrix X 
is small. Let F — KC := T so that A + BKC = A + BF — BT. Then using the 
formula in [12] we have 

|A i(A + BKC) - A j(A + BF ) | < ||BT|| 2 *P0 

< ||fl|| a ||T|| 2 *(X) (3.9) 

<||B|| 2 ||F-.ffCM(A') 

which shows that control over the eigenvalue locations of A + BKC can be obtained 
only if both \\F— KC\\ and k(X) are kept small. One way of doing this is to minimize 

J = oii(X) + o 2 ||F-KC|||. 

= ai ^)Al+a 2 Tt!ux{(F-KC) T (F-KC)} (310) 

by letting A (A + BF) range over the region 0 C C~. Similarly, by letting 
A + FdC = A + BKdC a dual problem can be formulated as 

Jd = fix + /? 2 Trac e{(F D - BI< D ) T (F D - BK D )} (3.11) 

GminyA-D) 

The idea of improving the conditioning of the eigenstructure and of minimizing 
the norm of F — KC was first introduced in Keel and Bhattacharyya [13] [14] . Here 
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an improved version of this algorithm is presented. In particular we convert the 
constrained optimization problem to an unconstrained problem and extend the class 
of regions 0, C C~ to more general and useful regions. These details are given next. 

4. STABILIZATION ALGORITHM 

In the Sylvester equation approach described in Section 2, 

AX - XA = —BG (4.1) 

FX = G (4.2) 

and let A(A) C C C~. Under the assumption A(A) D A(A) = 0 and (A,B) 

controllable, ( G , A) observable, the unique solution X will ‘almost surely’ be non- 

singular by deSouza and Bhattacharyya [7] and then A(A + BF) = A(A) with 
F = GX~ l . By letting A(A) range over Q this algorithm generates the family of 
F(Cl), by letting G be a free parameter run through all possible values this formula 
generates the family F(fl) defined in (3.8). 

If A is a complex diagonal matrix in (4.1), it is clear that X in (4.1) is the 
corresponding complex eigenvector matrix. However we want to treat these matrices 
as real for computational convenience. The following Lemma 4.5 shows that A can 
be taken as a real matrix without loss of generality. Before we state Lemma 4.3 it 
is necessary to introduce some facts. 


Definition 4.1 
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A real square matrix A is called a pseudo diagonal matrix if it is of the form 

<*1 01 
-01 oil 

02 02 
~02 Ol 2 

Oi 3 

with o,-, 0i real. 

Definition 4.2 

A complex square matrix is called normal if A* A = AA* . 

Lemma 4.3 [15] 

A complex square matrix is unitary similar to a diagonal complex matrix if 
and only if it is normal. 

Lemma 4.4 




Any real pseudo diagonal matrix is normal. 
Proof 

Taking the i th block from (4.3) such as 





(4.4) 

(4.5) 


A = Diag(Ai A 2 A n ) 


(4.6) 
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AA* = Diag( A! A; 

A 2 A\ • • • A n A* ) 

(4.7) 

A* A = Diag( A*Ai 

A* 2 A 2 ••• A* n A n ) 

(4.8) 


Since AA* = A* A, the statement is true. 


Lemma 4.5 

Let (A + BF)X = XA and (A + BF)Y = YA where 

1. A, B, A, X and F are real matrices with appropriate dimensions. 

2. A is real pseudo diagonal, .4 is complex diagonal, and 

3. X and Y are nonsingular. Then, 

k(X) = k(Y) (4.9) 


Proof 

From Lemma 4.1 and 4.2, A is known to be normal and unitary similar to the 
complex diagonal matrix A. Thus 


A = uAu*. 


(4.10) 


Write 


so that 


and 


Now, 


(A + BF)X = XA = XUAU* 

(A + BF)XU = XU A 
XU = Y. 

YY' = XUU’X' = XX" = XX T . o 


(4.11) 

(4.12) 

(4.13) 

(4.14) 
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From this Lemma, minimizing cr max (X) / cr m in{X) in (3.10) is equivalent to 
minimizing cr mai (F)/o' m , n (} /r ). Therefore we can henceforth take A as a real pseudo 
diagonal matrix without loss of generality. In fact the condition numbers of X and Y 
are equal, i.e. k(X) = k(XU ) = k(Y). In order to use a gradient based algorithm 
the closed form expression of the gradient of the performance index (3.10) with 
respect to the variables G , K and the variable elements of A denoted a,- is evaluated. 
The details of this derivation are given in Appendix A. 

Theorem 4.6 


Given the performance index J in (3.10), and constraints (4.1) and (4.2), the 
gradients of J with respect to the independent variables G, K, and A are as follows: 

(a) 

= 2{a 2 (F - KC)X~ t + B t U t } (4.15) 

where U satisfies 


~ Ot j 

AU U A — j / -tr\ { < 7min(X')V a U a ^roai(-Y)l , iUj } 

-2a 2 X~ 1 (F T ~(KC) t )F 


(4.16) 


where v a and u a are right and left singular vectors corresponding to cr max (A?) 
and Vi and u; are for cr m i n (X ), respectively. 

(b) Let a, denote a variable element of A: 


dJ 

da,i 


-Tra.ce{UX^} 


(4.17) 


(c) 


where U satisfies (4.16) 


6 J_ 

dK 


-2 a,(F - I<C)C T 


(4.18) 
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Equations (4.15) — (4.18) are used to devise a gradient algorithm that iterates 
on the free parameters G, K and the entries of A to reduce J. At each iteration 
of the algorithm we get Ai,Fj and A^. Since cr(A { ) C Cl we have a(A + BFi) C Cl 
for each i. However a(A + BK{C ) may or may not be in Cl for each i, and the 
algorithm is designed to make a(A + BK{C) close to cr(A;) = a (A + BF t ) after 
some iterations. 

The following structure of the closed loop eigenvalue matrix A ensures stability 
without constraints during the iterations. 

a\ — a\ 

Note that a.i in the matrix A are the only nonzero parameters and furthermore 
the stability requirement, <r(A) C C~, can be automatically satisfied without 
constraints, for all real values of 5;. 

We can also parameterize A in such a way that the desired closed loop 
eigenvalue locations axe automatically confined to some useful region Cl as in Figures 
4.1 and 4.2. 

In choosing A, a maximal number of 2x2 blocks axe included in the initial choice. 
As the algorithm evolves some of the off diagonal terms may become very small. 
At that point we start to vary the corresponding diagonal terms independently. In 
the damping ration region described in Figire 4.2. 6 is also a free parameter. 
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Marginal Stability Region 

For this case we can simply modify the matrix A to 

/-(af+7) a 2 

-a 2 ~(a\ + 7) 

-(« 3 + 7) «4 

7 -04 -(5§ + 7) 


V 


\ 




with a,- as the real variable parameters and 7 is fixed. The eigenvalues of A are all 
to the left of the line Re(s ) = —7. 


Damping Ratio Region 


/ — (5j + 7) (a| + 7)tan^»sin^i 

[ — (Sj + 7)tan^>sin^i — (af + 7) 


V 


\ 




Now we discuss what happens when the proposed algorithm fails to find a 
stabilizing controller of order i. In this case, we increase the controller order to 
t + 1. It is then necessary to have a way to select the initial values of Go, Ao and 
Kq for the controller of order i + 1 to ensure that the performance index J keeps 
decreasing. The following theorem shows the way to select initial variables so that 
J always decreases with increasing controller order. 

Theorem 4.7 

Let J* be the optimal performance index with optimal variables G* , A* and 
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K* where 


and X* and F* satisfy 


j. = + | jF . _ K-cfr 


AX * - X*A* = -BG* 


f* = g*(x* r 1 


( 4 . 19 ) 


Then for the extended system 


A 0 


B e = 


B 0 
0 Ii 


C ( 


_(C 0 

V° h. 


( 4 . 20 ) 


e V° IiJ 

the value of its performance index J t is equal to J* if the set of initial variables are 


Ge — 


(G 0 \ _(K _0 \ 

VO X 3 Ai) Ae_ V0 XzAiX; 1 ) 


0 X 3 A { 

where 4, is an arbitrary pseudo diagonal matrix of a extended matrix 


( 4 . 21 ) 


-(■: D 


and 




X 3 = 


\ 


v 2 


V Oi) 

for <Ti > tr 2 > • ■ • > > 0 with cr» > <T mtn (X*) and 0 \ < cr max (X*). 


Proof 


Let the optimal values of J* be obtained by G* and A*, then the extended 
system becomes 

(a o \(X' Xi\_(x* XA(A* o\ 

Vo 0 J\X 2 . X 3 ) \X 2 X 3 J \0 AiJ 



I 



I 
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__(B 0 \(G' GA 

\ 0 Ii J \ G 2 G 3 J 

Be G c 

( AX*-X*A * AXi - XiAA ( BG* BGA 

V -X 2 A -X 3 Ai J \ G 2 g 3 J 


if we pick G\ = 0 and G 2 — 0, then X x = 0 and X 2 - 0 and X 3 A, 
choose 


I*' 


\ 


<?i 


) 


G 3 . Here we 


(4.23) 


for < 7 i > <r 2 > • • • > < 7 ,- > 0 with a ,• > cr m i n (X*) and a x < cr max (X*). Such a X 3 is 
guaranteed by the choice of G 3 = X 3 A{ and 


G min 


(X*) = <r m ,n(X e ) 


(4.24) 


G max 


(X*) 


.(*<) 


Therefore, 


Gjnax (jn G max (JQ 

&min (X*) 

Grain (Xe) 


Now consider the term ||.F — KC\\ 2 F . Since 


we have 


where 



(4.25) 


(4.26) 


(4.27) 



I 


Now 


G 0 \ ( X- 1 0 

e - - l o X 3 Ai Mo X- 1 


F e = G e X7 l = 


GX 


T — 1 


and let 


then 


AT = 


0 X 3 A i X7 x 


K I\\ 


K 2 k 3 


F _ Kr _ (gx-'-kc -K x 
Fe 1 eCe - V -K 2 c X 3 AiX 3- 1 - k 3 

Here we choose I\\ = 0 and K 2 = 0. Also we can choose 


K 3 = X 3 AiX7 l 

because X 3 and A t are well defined. With such a K we have 

GX~ x -KC 0\ 


F e - K e C e = 


0 


0 


Thus, 


\F - KC\\ 2 F = \\F e - K e C e \\ 2 F 


Therefore, we conclude 


<T ” , " ( y* ) + IIJ” - K'cnl = + II .F. - K'C. fr 

&min\X ) (J m in\X e ) 


with choices of 


G* 0 


K 


0 


and K e - 1 Q x 3 AiXf 1 


0 X 3 Ai 

with X 3 as in (4.23). This concludes the proof. o 
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(4.28) 

(4.29) 

(4.30) 

(4.31) 

(4.32) 

(4.33) 

(4.34) 

(4.35) 


This theorem is useful for finding a low order stabilizing controller because it 
shows how by sequentially increasing the order of the controller, J can be guaranteed 
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to decrease. Since a small enough value of each term of J confines the spectrum of 
-4 + BKC to Q (in accordance with (3.9)) the algorithm eventually stabilizes the 
system by sequentially increasing the order of controllers. 

5. EXAMPLES 

The algorithm developed in the last section is applied to several examples here. 
The gradient calculations of Theorem 4.7 are used along with the Harwell subroutine 
package. 

Example 1 

The first example is a simplified model of the NASA F -8 Digital Fly-By- 
Wire(DFBW) airplane[16] and its dynamic equation of lateral directional is as 
follows. 

/ p\ / —2.6 0.25 —38. 0 \ /P\ / 17. 7. \ 

d_ l r j _ [ -0.075 -0.27 4.4 0 ] [ r ] [ 0.82 -3.2 / 6 a \ 

dt P ~ 0.078 -0.99 -0.23 0.052 /? + 0 0.046 I \S r ) 

\<t>) \ 1.0 0.078 0 0 / \<t>) \0 0 / 

(r\ _( 0 1 0 0 \ [ r | 

\<t>) - v° o o i ) yj 

The given design specifications [16] are that the closed loop poles must be the left 
of the line s = —0.2 i.e. 7 = 0.2, and the damping factor is > 0.7, i.e. <j> = j in 
Figure 4.2. total equilibrium velocity Vo = 620 ft/s(Mach = 0.6) and equilibrium 
angle For the optimization problem initial values are chosen to be 
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/I 1.5 0.5 — 2\ 

0 \5 1 -0.25 0.5 ) 

A '° = (o o) 

After 41 gradient iterations minimizing J in (3.10) the following 0 th order 

stabilizing compensator is obtained. 

/ 4.60357 — 1.75629 \ 

V 5.21515 —1.85922 J 

Note that the order of pole placement compensators (both Brasch - Pearson and 
Kimura) is 1. The corresponding data is shown in Tables 1.1, 1.2 and Figure 5.1. 
For comparison, the same problem was run without including the condition number 
term in J (i.e. ai = 0 in (3.10)). It is seen from the corresponding data, shown in 
Table 1.3, 1.4 and Figure 5.2 that the condition number increases significantly, and 
although stabilization is achieved the closed loop eigenvalues fail to be in Q. 


TABLE 1.1 
Eigenvalues 

(<*i = 1, a 2 = 1, <f> = f , C > 0.7, 7 = -0.1) 


Aq 

A 0 + BqKqC 

Ao + B 0 Fq 

Aq + BoKqCo 

—2.39 ± j'0.00 
+0.00 ± jO.OO 
-0.34 ± >2.62 

-2.39 ± >0.00 
+0.00 ±>0.00 
-0.34 ±>2.62 

-9.45 ± >3.70 
-0.34 ± >0.29 

-7.58 ± >4.96 
-0.42 ±>0.33 


TABLE 1.2 
Objective Function. 



J 

mmamwm 

k(X 0 ) 

Initial 

Optimal 

155.9021 

47.03439 

61.3301 

0.06839 

94.572 

46.966 
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TABLE 1.3 
Eigenvalues 

(<*i=0, <*2 = l,^=f,C>0.7, 7 = 0.1) 


^•0 

A 0 +B 0 K°C 0 

A 0 + BqFq 

A 0 + B 0 K*Co 

-2.39 ± iO.OO 
+0.00 + >0.00 
—0.34 ± j2.62 

-2.39 ± j'0.00 
+0.00 ± jO.OO 
—0.34 ±j2.62 

-2.39 ± jO.Ol 
-2.42 ± j'0.32 

— 1.44 ±j 2.54 
-3.44 ± >0.00 
-1.15 ± j'0.00 


TABLE 1.4 
Objective Function 



J 

\\F 0 - K 0 Co\\ 2 f 

k(X 0 ) 

Initial 

Optimal 

155.9021 

233089.02 

61.3301 

0.01530 

94.572 

233089 


Example 2 


Consider the symmetric vibration model of the standard Draper/RPL satellite 


shown in Figure 5.3. The dynamic equations, taken from [17] are: 







20 


where 


A = 


/° 

0 

0 

0 

0 


0 

0 

0 

14.8732 

-146.702 


VO -41.8468 

/ 


0 

0 

0 

32.8086 

-7476.64 

-2699.36 


1 0 0 \ 
0 
1 
0 
0 

0/ 


0 1 
0 0 
0 0 
0 0 
0 0 


B = 


0 

0 

0 

0.04168 


0 

0 

0 

0.23623 


\ 


10.38611 -25.647 

V 3.725120 -9.1629/ 


C = 


0 10 0 
0 0 10 


1 0> ) 

0 \) 


From the design specifications in [17], it follows that the closed loop system 
must have poles to the left of s = —0.5. For the minimization of J the initial values 
are chosen to be 


! — 0.2 2 
-2 - 0.2 


A 0 = 


\ 


-1 10 

-10 -1 


-0.5 


„ / 1.125 1. 

G ° “ \ -1 2 . 


1.5 

-0.5 

3.5 

-1 -0.5 / 

1.5 2 \ 

2.5 

1.6 

4 

0.5 -1 J 

K 0 = 

/ 0 0 
lo 0 

) 



After 67 iterations, the following 0 th order stabilizing controller is obtained: 


-90.97491 20.62868 \ 

A “ V -197.646 5.326668 J 

Note that the order of pole placement compensators (both Brasch - Pearson[3] and 
Kimua[4]) is 3. Tables 2.1, 2.2 and Figure 5.4 display the performance indices and 
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the corresponding eigenvalue locations. For the purpose of comparison, the problem 
was also run with the condition number term left out of the performance index (i.e. 
a i = 0). In this case the algorithm fails to stabilize the system as shown in Table 
2.3, 2.4 and Figure 5.5. This example illustrates that both terms of the performance 
index need to be considered in the stabilization procedure. 


TABLE 2.1 

Eigenvalues for Example 2 
(<*i = 1, a 2 = 1, T = 0.5 ) 


Ao 

A 0 + B 0 K°C 

Ao + BoFq 

Ao + BqKq Co 


+0.00 ± >53.1 

+0.00 ±>53.1 

-2.89 ±>36.7 

-3.58 ±>30.7 


+0.00 ±>5.43 

+0.00 ± >5.43 

-2.18 ±>0.30 

-2.41 ±>0.67 


+0.00 ±>0.00 
+0.00 ±>0.00 

+0.00 ± >0.00 
+0.00 ± >0.00 

-0.88 ± >5.82 

-1.45 ±>6.08 



TABLE 2.2 
Performance Indices 



J 

mwamm 

k(Xo) 

Initial 

Optimal 

11965915 

587.2232 

11965506 

45.63520 

409.2925 

541.5880 


TABLE 2.3 

Eigenvalues for Example 2 
(<*o = 0, a 2 = 1, 7 = 0.5) 


^4.0 

Ao + BqKqCo 

Ao + BqFq 

Ao + BqKqCo 

+0.00 ±>53.1 
+0.00 ±>5.43 
+0.00 ±>0.00 
+0.00 ±>0.00 

+0.00 ±>53.1 
+0.00 ± >5.43 
+0.00 ± >0.00 
+0.00 ±>0.00 

-0.63 ±>0.05 
-0.66 ± >3.41 
-0.59 ±>7.86 

+178. ±>0.00 
-2.57 ±>6.19 
+1.61 ±>0.00 
+0.83 ±>1.07 
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TABLE 2.4 
Performance Indices 



J 

\\Fo-K 0 Co\\ 2 f 

k(X 0 ) 

Initial 

Optimal 

11965915 

383859.66 

11965506 

490.5633 

409.2925 

383369.1 


6. CONCLUDING REMARKS 

The results given here are algorithmic in nature and can be improved upon by 
developing constructive necessary and sufficient conditions for stabilizability with a 
fixed order controller. This in turn will require effective ways of characterizing 
the Hurwitz region. These problems are difficult and have received very little 
attention in the literature. Finally we mention that the algorithm guarantees neither 
a “global” minimum nor does it always find a stabilizing controller of a prescribed 
order whenever one exists. The existence of stabilizing controllers of a fixed order 
is still our unsolved problem. 
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APPENDIX 


Proof of Theorem 3.4 

(a) 


J = a i + a 2 Trace{(F - KC) T (F - KC )} 


&min (X) 


Let 


J i : = 


=(-*■) 


.(*) 


= Trace { — -} 


and 


Trace { ^ (<r m i n (X)A<T mar (A )) cr maiA<T mtn (X)} 

°m.nl A J 


Note that 


A<7 max (X) — w a AXr a 


A<x m i n (X) — Mj AXuj 


where u,- and u, are left and right singular vectors corresponding to a 
and v a are for (T max . Thus, 


AJ] = 




Trace{ cr Tn i„(X)viu i }A.X 


Now 


J 2 : = Trace{(F - I<C) T (F - KC)} 

= Trac e{F T F - ( KC) T F - F T (KC) + (KC) T (KC)} 

= Trace(E T F) - 2Trace{(KC) T F} + Trace{(AC) T (/YC)} 


(-4.1) 

(-4.2) 

(A. 3) 

(- 4 . 4 ) 

(-4.5) 

and t; a 

M-6) 

(A.7) 


A J 2 = 2Trace(F T AF) - 2Trace{(/\ C) T AF} 
= 2Trace{[F r - (KC) T ]AF) 


Now we have 


AJ = 






Trace{cr m ,„(X)ujU i (T max {X')v a u r £} AX 

+ 2a 2 Trace{(F r - ( I<C) T )AF }. 


From F = (?JY 1 , the gradient of F with respect to G is given directly as 


AF = AGX' 1 + GA(X -1 ) 

= AGJY" 1 -GX- 1 AXX~ l 
= A GX- 1 - FAXX~ l 
= {AG- FAX)X~\ 


Substituting (.4.10) into (A. 9) we have 


A J = 2a 2 Trace{(F r - {KC) T )AGX~ 1 } 

C^i 

4” Trace{ ^ y^r(o' m i n {X')viU^ &max{X')v a u a ) 

- 2 a 2 X~ 1 (F T - {KC) t )F}AX 


Using [14] 

AAX - AX A = -BAG 

and 

. n n 

AA = 5Z Tij4'~ 1 .BAG4 '~ 1 . 


(.4.8) 


(4.9) 


(4.10) 


(4.11) 


(4.12) 


(4.13) 
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Substituting (.4.13) into the second term of (4.11) we have 


A J : = 2o 2 Trace{A r-1 (F r - (KC) 1 )AG} 

n n 

+ Trace{J] J] jijA^ 

i — i i 


( — -^ I (iKt i :)-2« 2 r , (F J -(AXT)^) 

a min l < A ) 


A i ~ 1 BAG} 

= 2a 2 Trace{X- 1 (F r - (7i'C) T )AG} 

n n 

+ Trace { ^ ^ ■jijA 1 ’ -1 X/A 1-1 BAG} 

1=1 j = l 

S ^ 

u 

= Trace{[2a 2 X -1 (F r - ( KC) T ) + BU))AG 
From (4.12) and (4.13) it follows that U is the unique solution of 

AU ~ U A = X f . 


(4.14) 


(4.15) 


Therefore 


where U satisfies 


= 2{a 2 (F - KC)X~ t + B T 77 r } 


4C7 - £74 = 


(4.16) 


<*i 


{<T mtn (X>,u/' - <r m a X (X)v a uZ} _ 2a 2 X" 1 (F r - ( KC) T )F (4.17) 


° 2 min(X) 


(b) 


Now we evaluate the gradients of (3.10) with respect to the variable elements 


of 4. Recall the equation (4.9) 
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AJ = 




Trace{cr m i n (X)v a u a c7 maj; (A)vj'Uj }AX 


1 77 min 2 (X) 

+ 2a 2 Trace{( J F r - ( I<C) T )AF }. 

From F — GX ~ 1 , we compute( G is fixed) 

AF = -GX~ 1 AXX ~ 1 

= -FAXX~\ 


(A18) 


(A19) 


Substituting AF into (A. 18) 


AJ — Trace 


{ 2 — T^ri^miniX^iuJ - (T max (X :)v a ul) - 2 a 2 X *(F - (KC) t )F} AX 

a minV A ) 


T, 


= Trace {X/ A A} 


Since 


AAX - AX A = XAA 
AX = ££>, ‘ A ‘-'(- XAA ) Ai ~' 

i= 1 J=1 


(A.20) 

(A.21) 
(A. 22) 


Substituting (A. 22) into (A.20) 


AJ = Trace} ^ ^ 7 ij x /A ,_1 (— XAA)A J_1 } 

t=i j=i 
n n 

= -Trace} £ ^^A^X/A'" 1 XAA} 


(A. 23) 


«=i >=i 


1 / 


It is clear that U is the unique solution of 


AU -U A = X f 
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as in (A 14). 


Therefore, 


A J = -Tiace{UXAA} 


dJ _ (TT „dA. 
^ = -Trace{tOC— } 


As an example the following calculation is considered. Let 


U = 


U 11 U 12 

U 21 u 22 


X = 


X 11 X i2 
X 21 x 22 


A=l a j °2 

0 a\ 


Then 


— -2Tra ce I ( Un Ul2 ^ 
da, 2T IU21 «2 2 ) 

= 4 a,(u„x n + v-i2 x 2i) 

2 =2Trace(f““ “ 12 ) 
a 2 ^ \ U 21 U 22 J 

= 4 a 2 (u 2 ixi 2 +U22Z22) 


dJ 

da 


x \\ x 12 
X 21 x 22 


x l\ X, 2 
X 21 x 22 


2a 1 0 

0 0 


0 0 
0 2a 2 


or 



_ 4 / di(una:ii + rti 2 :r 2 i) 
\a 2 (u2lXi2 + U 2 2 x 22 ) 


(c) 


Finally the gradient of J with respect to K is easily derived. 

A J = -2a 2 Trace{CF T AJv - C(KC) T AK} 

= — 2a 2 Trac e{(CF T - C(KC) T )AK] 


(A. 24) 
(A. 25) 

(A. 26) 

(A. 27) 
(A. 28) 

(A. 29) 
(A. 30) 


Thus, 


dJ 

dK 


-2a 2 [e - KC]C T 


(A31) 
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